{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# Reference: https://jupyterbook.org/interactive/hiding.html\n",
    "# Use {hide, remove}-{input, output, cell} tags to hiding content\n",
    "\n",
    "import sys\n",
    "import os\n",
    "if not any(path.endswith('textbook') for path in sys.path):\n",
    "    sys.path.append(os.path.abspath('../../..'))\n",
    "from textbook_utils import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "csv_file = 'data/Full24hrdataset.csv'\n",
    "usecols = ['Date', 'ID', 'region', 'PM25FM', 'PM25cf1', 'RH']\n",
    "full = (pd.read_csv(csv_file, usecols=usecols, parse_dates=['Date'])\n",
    "        .dropna())\n",
    "full.columns = ['date', 'id', 'region', 'pm25aqs', 'pm25pa', 'rh']\n",
    "\n",
    "bad_dates = ['2019-08-21', '2019-08-22', '2019-09-24']\n",
    "GA = full.loc[(full['id'] == 'GA1') & (~full['date'].isin(bad_dates)) , :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "from scipy.stats import randint\n",
    "\n",
    "X = GA[['pm25aqs', 'rh']]\n",
    "y = GA['pm25pa']\n",
    "\n",
    "theta2_hat = LinearRegression().fit(X, y).coef_[1]\n",
    "\n",
    "def boot_stat(X, y):\n",
    "    n = len(X)\n",
    "    bootstrap_indexes = randint.rvs(low=0, high=(n - 1), size=n)\n",
    "    theta2 = (\n",
    "        LinearRegression()\n",
    "        .fit(X.iloc[bootstrap_indexes, :], y.iloc[bootstrap_indexes])\n",
    "        .coef_[1]\n",
    "    )\n",
    "    return theta2\n",
    "\n",
    "np.random.seed(42)\n",
    "boot_theta_hat = np.array([boot_stat(X, y) for _ in range(10_000)])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basics of Confidence Intervals\n",
    "\n",
    "We have seen that modeling leads to estimates, such as the typical time that a bus is late ({numref}`Chapter %s <ch:modeling>`), a humidity adjustment to an air quality measurement ({numref}`Chapter %s <ch:linear>`), and an estimate of vaccine efficacy ({numref}`Chapter %s <ch:data_scope>`).  These examples are point estimates for unknown values, called *parameters*: the median lateness of the bus is 0.74 minutes; the humidity adjustment to air quality is 0.21 PM2.5 per humidity percentage point; and the ratio of COVID infection rates in  vaccine efficacy is 0.67. However, a different sample would have produced a different estimate. Simply providing a point estimate doesn't give a sense of the estimate's precision. Alternatively, an interval estimate can reflect the estimate's accuracy. These intervals typically take one of two forms:\n",
    "\n",
    "1. A *bootstrap confidence interval* created from the percentiles of the bootstrap sampling distribution\n",
    "1. A *normal confidence interval* constructed using the standard error (SE) of the sampling distribution and additional assumptions about the distribution having the shape of a normal curve\n",
    "\n",
    "We describe these two types of intervals and then give an example."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that the sampling distribution (see {numref}`Figure %s <triptych>`) is a probability distribution that reflects the chance of observing different values of $\\hat{\\theta}$. Confidence intervals are constructed from the spread of the sampling distribution of $\\hat{\\theta}$, so the endpoints of the interval are random because they are based on $\\hat{\\theta}$. These intervals are designed so that 95% of the time the interval covers $\\theta^*$.   "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As its name suggests, the percentile-based bootstrap confidence interval is created from the percentiles of the bootstrap sampling distribution. Specifically, we compute the quantiles of the sampling distribution of $\\hat{\\theta}_B$, where $\\hat{\\theta}_B$ is the bootstrapped statistic. \n",
    "For a 95th percentile interval, we identify the 2.5 and 97.5 quantiles, called $q_{2.5,B}$ and $q_{97.5,B}$, respectively, where 95% of the time the bootstrapped statistic is in the interval:\n",
    "\n",
    "$$ \n",
    "q_{2.5,B} \\leq \\hat{\\theta}_B~ \\leq ~ q_{97.5,B}\n",
    "$$\n",
    "\n",
    "This bootstrap percentile confidence interval is considered a quick-and-dirty interval. There are many alternatives that adjust for bias, take into consideration the shape of the distribution, and are better suited for small samples. \n",
    "\n",
    "The percentile confidence interval does not rely on the sampling distribution having a particular shape or the center of the distribution being $\\theta^*$. In contrast, the normal confidence interval often doesn't require bootstrapping to compute, but it does make additional assumptions about the shape of the sampling distribution of $\\hat{\\theta}$. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the normal confidence interval when the sampling distribution is well approximated by a normal curve. For a normal probability distribution, with center $\\mu$ and spread $\\sigma$, there is a 95% chance that a random value from this distribution is in the interval $\\mu ~\\pm ~ 1.96 \\sigma$. Since the center of the sampling distribution is typically $\\theta^*$, the chance is 95% that for a randomly generated $\\hat{\\theta}$: \n",
    "\n",
    "$$|\\hat{\\theta} -\\theta^*| \\leq 1.96 SE(\\hat{\\theta})$$\n",
    "\n",
    "where $SE(\\hat{\\theta})$ is the spread of the sampling distribution of $\\hat{\\theta}$. We use this inequality to make a 95% confidence interval for $\\theta^*$:\n",
    "\n",
    "$$ [ \\hat{\\theta} ~-~ 1.96 SE(\\hat{\\theta}),~~~  \\hat{\\theta} ~ +~ 1.96 SE(\\hat{\\theta})]$$\n",
    "\n",
    "Confidence intervals of other sizes can be formed with different multiples of $SE(\\hat{\\theta})$, all based on the normal curve. For example, a 99% confidence interval is $\\pm 2.58 SE$, and a one-sided upper 95% confidence interval is $[ \\hat{\\theta} ~-~ 1.64 SE(\\hat{\\theta}),~~ \\infty]$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{note}\n",
    "\n",
    "The SD of a parameter estimate is often called the *standard error*, or SE, to distinguish it from the SD of a sample, population, or one draw from an urn. In this book, we don't differentiate between them. We call them SDs. \n",
    "\n",
    ":::"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We provide an example of each type of interval next."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Earlier in this chapter we tested the hypothesis that the coefficient for humidity in a linear model for air quality is 0. The fitted coefficient for these data was $0.21$. Since the null model did not completely specify the data generation mechanism, we resorted to bootstrapping. That is, we used the data as the population, took a sample of 11,226 records with replacement from the bootstrap population, and fitted the model to find the bootstrap sample coefficient for humidity. Our simulation repeated this process 10,000 times to get an approximate bootstrap sampling distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the percentiles of this bootstrap sampling distribution to create a 99% confidence interval for $\\theta^*$. To do this, we find the quantiles, $q_{0.5}$ and $q_{99.5}$, of the bootstrap sampling distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lower 0.05th percentile: 0.099\n",
      "Upper 99.5th percentile: 0.260\n"
     ]
    }
   ],
   "source": [
    "q_995 = np.percentile(boot_theta_hat, 99.5, method='lower')\n",
    "q_005 = np.percentile(boot_theta_hat, 0.05, method='lower')\n",
    "\n",
    "print(f\"Lower 0.05th percentile: {q_005:.3f}\")\n",
    "print(f\"Upper 99.5th percentile: {q_995:.3f}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, since the histogram of the sampling distribution looks roughly normal in shape, we can create a 99% confidence interval based on the normal distribution. First, we find the standard error of $ \\hat{\\theta} $, which is just the standard deviation of the sampling distribution of $\\hat{\\theta}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.02653498609330345"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "standard_error = np.std(boot_theta_hat)\n",
    "standard_error"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, a 99% confidence interval for $\\theta^*$ is $2.58$ SEs away from the observed $\\hat{\\theta}$ in either direction:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lower 0.05th endpoint: 0.138\n",
      "Upper 99.5th endpoint: 0.275\n"
     ]
    }
   ],
   "source": [
    "print(f\"Lower 0.05th endpoint: {theta2_hat - (2.58 * standard_error):.3f}\")\n",
    "print(f\"Upper 99.5th endpoint: {theta2_hat + (2.58 * standard_error):.3f}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These two intervals (bootstrap percentile and normal) are close but clearly not identical. We might expect this given the slight asymmetry in the bootstrapped sampling distribution. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are other versions of the normal-based confidence interval that reflect the variability in estimating the standard error of the sampling distribution using the SD of the data. And there are still other confidence intervals for statistics that are percentiles, rather than averages. (Also note that for permutation tests, the bootstrap tends not to be as accurate as normal approximations.) "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{note}\n",
    "\n",
    "Confidence intervals can be easily misinterpreted as the chance that the parameter $\\theta^*$ is in the interval. However, the confidence interval is created from one realization of the sampling distribution. The sampling distribution gives us a different probability statement; 95% of the time, an interval constructed in this way will contain $\\theta^*$. Unfortunately, we don't know whether this particular time is one of those that happens 95 times in 100 or not. That is why the term _confidence_ is used rather than _probability_ or _chance_, and we say that we are 95% confident that the parameter is in our interval. \n",
    "\n",
    ":::"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Confidence intervals and hypothesis tests are related in the following way. If, say, a 95% confidence interval contains the hypothesized value $\\theta^*$, then the $p$-value for the test is less than 5%. That is, we can invert a confidence interval to create a hypothesis test. We used this technique in the previous section when we carried out the test that the coefficient for humidity in the air quality model is 0. In this section, we have created a 99% confidence interval for the coefficient (based on the bootstrap percentiles), and since 0 does not belong to the interval, the $p$-value is less than 1% and statistical logic would lead us to conclude that the coefficient is not 0. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another kind of interval estimate is the prediction interval. Prediction intervals focus on the variation in observations rather than the variation in an estimator. We explore these next. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
